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psj , Abstract 

Recovery algorithms play a key role in compressive sampling (CS). Most of current CS recovery 
I algorithms are originally designed for one-dimensional (ID) signal, while many practical signals are two- 

. dimensional (2D). By utilizing 2D separable sampling, 2D signal recovery problem can be converted into 

ID signal recovery problem so that ordinary ID recovery algorithms, e.g. orthogonal matching pursuit 
(OMP), can be applied directly. However, even with 2D separable sampling, the memory usage and 
Q I complexity at the decoder is still high. This paper develops a novel recovery algorithm called 2D-OMP, 

which is an extension of ID-OMP. In the 2D-0MP, each atom in the dictionary is a matrix. At each 
^ . iteration, the decoder projects the sample matrix onto 2D atoms to select the best matched atom, and 

' then renews the weights for all the already selected atoms via the least squares. We show that 2D-OMP 

is in fact equivalent to ID-OMP, but it reduces recovery complexity and memory usage significantly. 
What's more important, by utilizing the same methodology used in this paper, one can even obtain higher 
dimensional OMP (say 3D-OMP, etc.) with ease. 
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I. Introduction 

Let a; G be a one-dimensional (ID) signal and 'J' € M"^" an orthonormal transform matrix, where 
M is the set of real numbers. If a; = ^'z and there are only k <^n spikes (nonzero entries) in z, we say 
that X is /c-sparse in ^ domain. We sample by $ € M™^" {k < m < n) to get y = = Az € 
where A = If $ obeys the order-Zc restricted isometry property (RIP) and has low coherence with 
then z (and in turn x) can be effectively recovered from y HI, ||2l, IJl. Many algorithms have 
been proposed to recover x from its random sample y, e.g. linear programming (LP) H and orthogonal 
matching pursuit (OMP) Q. For a detailed overview on recovery algorithms, please refer to IS. 

In practice, many signals, e.g. image, video, etc, are two-dimensional (2D). A straightforward imple- 
mentation of 2D compressive sampling (CS) is to stretch 2D matrices into ID vectors. However, such 
direct stretching increases exponentially the complexity and memory usage at both encoder and decoder. 
An alternative to ID stretching is to sample rows and columns of 2D signals independently by using 
separable operators |7 |. Through 2D separable sampling, encoding complexity is exponentially reduced. 
However, as the recovery problem is converted into a standard ID -minimization problem, decoding 
complexity is still very high. 

As a representative sparse signal recovery algorithm, the OMP achieves good performance with low 
complexity. The OMP is originally designed for ID signal recovery. To reduce the complexity of 2D 
signal recovery, this paper extends the ID-OMP to obtain the 2D-0MP. We prove that with 2D separable 
sampling, 2D-0MP is in fact equivalent to ID-OMP, so that both algorithms will output exactly the same 
results. However, the complexity and memory usage of 2D-0MP is much lower than that of ID-OMP. 
Thus, 2D-0MP can be used as an alternative to ID-OMP in 2D sparse signal recovery. 

This paper is arranged as follows. Section |ll] first briefly reviews the principles of 2D separable sampling 
and ID-OMP, and then makes a detailed analysis on the complexity of ID-OMP. In Section Hill we deduce 
the 2D-0MP algorithm, reveal the equivalence of 2D-0MP to ID-OMP, and compare the complexity and 
memory usage of 2D-0MP with that of ID-OMP. In Section |IVl simulation results are reported. Finally, 
Section |V] concludes this paper. 

II. ID Orthogonal Matching Pursuit with 2D Separable Sampling 

A. 2D Separable Sampling 

The principle of 2D separable sampling Q is as follows. Let X G M"^" be a 2D signal which is 
fc-sparse in ^ domain, i.e. X = ^fZ*^ and there are only k <^v? spikes in Z, where (•)^ denotes the 
transpose. For simplicity, we use the same operator $ to sample the rows and columns of X independently 
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to get Y = *X*'^ = AZA^ e M™x™. Let y G M™' be the ID stretched vector of Y and ^ G M"' the 
ID stretched vector of Z. It was proved that 



where (g) denotes the Kronecker product Q. It is easy to prove (g) = (g) 

hence Q = A(E) A. Now this is just a standard ID sparse signal recovery problem which can be attacked 
by the OMR 

B. ID Orthogonal Matching Pursuit 

Let Q, = (cji, ■ ■ ■ ,u>n^), where u-i G M™^ is the i-th column of fi. We call the dictionary and Ui 
an atom. The main idea of ID-OMP is to represent y as a weighted sum of as few atoms as possible. 

Algorithm [T] gives main steps of ID-OMP. To implement ID-OMP, we need two auxiliary variables. 
First, to avoid atom reselection, set i is defined to record the indices of those atoms that are allowed to 
be selected in the future (excluding those already selected atoms). Second, vector r G W^^ is defined to 
hold the residual after removing the selected atoms from y. Initially, r is set to y. Then at each iteration, 
the decoder picks from the dictionary the atom that best matches the residual and then renews the weights 
for all the already selected atoms via the least squares. 

C. Complexity Analysis 

We decompose the iteration of ID-OMP into the following steps and analyze its complexity step by 
step. 

1 ) Project: The projection of residual r onto atom uji> is (r, Ui') / Hg, where (•, •) denotes the inner 
product between two vectors and ||-||2 denotes the £2 -norm of a vector. Let p — (||u^x||r,,-- - , || tt'ns || 2) , 
then this step can be implemented by fl^r./p, where ./ denotes dot division. The complexity of this 
step is dominated by (n^ x m^) x (m^ x 1) matrix- vector multipUcation. Hence the complexity of this 
step is 0{'m?n?). 

2 ) Select Best Matched Unselected Atom: This step selects from unselected atoms the atom with the 
maximal absolute value of projection. As there are atoms and k <^n^, the complexity of this step is 
approximately O(n^), negligible compai^ed with the Project step. 

3) Renew Weights: Let fi^ = (o^j^, • • • ,u>i^), then Y^t'=i '^t'^if, = ^iU. According to linear algebra. 



where Q = QjQi G M*^* and g = i^Jy G M*. The complexity to calculate Q and g depends on t. For 
t < k <^ V?, the complexity of this step is negligible compared with the Project step. 



(1) 



argmin \\y — Qiu\\2 = Q g 



(2) 
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Algorithm 1: ID Orthogonal Mulching Pursuit with 2D Separable Sampling 
Input: 

• ft e R"*^^"^: sampling matrix 
. y G M"*': sample 

• k: sparsity level 
Output: 

» z E M"^: reconstruction of the ideal signal z 
Auxiliary Variables: 

• re W^^: residual 

• i: set of the indices of atoms that are allowed to be selected in the future 
Initialization: 

• r ^ y 

. i^{l,2,--- ,n^} 
for t 1 to A; do 

i ^ i\it; 

u ^ argmin ||y - 5]t'=i '"t'<^v II2' 
for t 1 to A; do 

Zit ^ Uu 



4) Update Residual: The complexity of this step depends on t. For t < A; ^ n^, the complexity of 
this step is negligible compared with the Project step. 

Based on the above analysis, we conclude that the complexity of ID-OMP is dominated by the Project 
step and its complexity is 0{m?'n?). 

III. 2D Orthogonal Matching Pursuit 

This section develops the 2D-OMP algorithm whose main idea is to represent 2D signal Y as a weighted 
sum of 2D atoms that are selected from an over-complete dictionary. We first redefine the concepts of 
atom, dictionary, and projection for 2D signals. Then we give the 2D-OMP algorithm. We reveal the 
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equivalence of 2D-OMP to ID-OMP and compare the complexity and memory usage of 2D-OMP with 
that of ID-OMP. 



A. Definition 

Let X G M"^" be a 2D signal that is /c-sparse in * domain and Y = ^X*^ = AZA^ e R"'^™. 
Let A = (ai, • • • , a„), where a, is the i-th column of A. We redefine dictionary, atom, and projection 
as follows. 

1 ) Dictionary and Atom: In the 2D-0MP, the dictionary contains atoms and each atom is an m x m 
matrix. Let Bj j be the (i, j)-th atom, then Bj is the outer product of and aj 



( 



B^j' — ciiCij 



auQij 



ana. 



\ 



(3) 



Now Y can be represented by the weighted sum of Bj j, i.e. 

n n 
1=1 3=1 

2) Projection: The projection of Y onto Bj j is 

(Y,B„-) 



IB 



(4) 



(5) 



where (Y,Bjj) = alYaj and ||Bjj ||2 is the Frobenius norm of Bjj, i.e. 



IB 



V ^ {ai,^aj'j)^ = ||ai||2 ||a 
\ i'=ij'=i 



J \\2 



(6) 



B. Algorithm Description 

Algorithm |2] gives main steps of 2D-0MP. To implement the 2D-0MP algorithm, we also need two 
auxiliary variables. First, to avoid atom reselection, set is defined to record the coordinates of those 
atoms that are allowed to be selected in the future (excluding those already selected atoms), where i for 
row indices and j for column indices. Second, R E K™x™ is defined to hold the residual after removing 
the selected atoms from Y. Initially, R is set to Y. Then at each iteration, the decoder first searches for 
the best matched atom in the dictionary and then renews the weights for all the already selected atoms 
via the least squares. 
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Alaorithm 2: 2D Orthogonal Malchinii Pursuit 



Input: 

• A G R"*^'*: sampling matrix 
, Y G R"»x"^: sample 

• k: sparsity level 
Output: 

• Z G M"^": reconstruction of the ideal signal Z 
Variable: 

. R G M"^x"*: residual 

• {i, j): set of the coordinates of atoms that are allowed to be selected in the future, i for row 
indices and j for column indices 

Initialization: 
. R^Y 

. (i j)^{(l,l),(l,2),-- - ,(n,n)} 
for f 1 to fc do 



(zt,jt)^arg max -^j- rH^; 

^ \ {itjt); 
w ^ arginin ||Y - '"t'B^,,j,, II2; 



u 



for t ^ 1 to A; do 



7 ) Least Squares: The weighted sum of t selected atoms constructs an approximation to Y. Let 

t 

K = Y-Y,Ui„j,Bi^„j,. (7) 

t'=i 

We model the problem as finding the optimal u = {ui^j^,--- ,Ui^jJ'^ that minimizes the Frobenius 
norm of R, which is in fact equivalent to the least squares problem. As ||R||2 = tr(RR^), where tr(-) 
is the trace of a matrix, the problem is equivalent to 

u = arg min tr(RR'^). (8) 
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Using dTJl and tr((-)"'") = tr(-), we have 

t 

tr(RRT) = tr(YYT)-2 J]ni,,,-,tr(YBT .J + 



t'=i 



t t 



i'=l s'=l 

||Y||2 + u'^Hw~2/^u, 



where 



I tr(B,,,,BT 



H 



tr(B,,,,,BT . ^ \ 



V tr(B.„,,BT^.J ... tr(B,„,,BT^.) j 



and 



/ = (tr(YB,^',. ),... ,tr(YBT .))T 



When tr(RR"'") takes the minimum, there must be 

atr(RRT) 



du 



2Hu - 2/ = 0. 



Hence 



« = H-i/. 

2j Calculation o/H <3?i(i /; It is easy to get 

Let (Bi^,j^,,Bj^,j^,) = {ai^,,ai^,) {aj^,,aj^,), then 

(Bji ji) j-j) 



H 



(Bii , Bjj 



Similarly, for /, because 



we have 



tr(YBT.) = tT{Ya,aJ) = ajYaj = (Y,B,j) , 



/ = ((Y,B,,,,Or- - ,(Y,B,,,,J)^ 



April 27, 2011 



8 



C. Equivalence of 2D-0MP to ID-OMP 



Obviously, the {n{i — 1) + j)-th atom of Q is 



auaj 



\ 



UJnii-l)+j 



(18) 



Compared with Q, it can be found that <^n(i-i)+j is just the ID stretched vector of Bj j. Hence, the 
Frobenius norm of Bj j will equal the ^2-norm of 

Then we prove that the projection of R onto Bj j equals the projection of r onto ijL)n{i-i)+j- Obviously, 



m m 



i'=lj'=l 




Hence 



(20) 



^nii-l)+j\\2 \\^i,j\\2 



It means: at each iteration of ID-OMP and 2D-0MP, the same atom will be selected. 

Finally, as r is the ID stretched vector of R, ||R||2 = ||'"||2- Hence, the least squares in ID-OMP and 
2D-OMP will output exactly the same results (in fact, it can be easily proved that Q = H and / = g). 

Based on the above analysis, we draw the conclusion that 2D-0MP is equivalent to ID-OMP. 

D. Complexity Analysis 

Below we analyze the complexity of 2D-0MP step by step. 

1) Project: Let P be an n x n matrix whose (i,j)-th element is ||Bjj||2. Then this step can be 
implemented by A^RA./P. The complexity of this step is dominated by (n x m) x (m x m) matrix- 
matrix multiplication and (n x m) x (m x n) matrix-matrix multipUcation. As m < n, the complexity 
of this step is 0{mv?). 

2) Select Best Matched Unselected Atom: Since there are n? atoms and k <C n^, the complexity of 
this step is approximately 0{v?), negligible compared with the Project step. 

3) Renew Weights: From ([141 ). it can be seen that the complexity to calculate (Bj^, j^, , Bj^, j^,) is 
0{m). Since (Y,Bjj) = ajYaj, the complexity to calculate (Y,Bj j) is 0{m?). The complexity to 
calculate H and / depends on t. For t < k <^ n?, the complexity of this step is negligible compared 
with the Project step. 
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Fig. 1. Speedup of 2D-0MP over ID-OMP. 



4) Update Residual: The complexity of this step depends on t. For t < k <^ -n?, the complexity of 
this step is negligible compared with the Project step. 

Based on the above analysis, we draw the conclusion that the complexity of 2D-0MP is 0{mv?), 
roughly 1/m of that of ID-OMP. 

E. Memory Usage Analysis 

For ID-OMP, an m? x v? matrix is needed to hold fi, so the memory usage is 0{m?v?), while for 
2D-0MP, fi is replaced by A, so the memory usage is reduced to 0{mn). 

IV. Experimental Results 

We have written both 2D-0MP and ID-OMP algorithms in MATLAB f8l. We present herein the 
results under three typical settings, i.e. (n, m) = (128, 16), (128, 24), and (128, 32). For each setting, we 
increase sparsity level k from 8 to 16. The transform matrix ^ is 128 x 128 2D discrete cosine transform 
(DCT) matrix. The sensing matrix $ is formed by sampling independent and identically-distributed (i.i.d.) 
entries from standard Gaussian distribution by using the randn function. 2D Sparse signal is obtained 
by using the sprandn function with density k/n'^. 

We run the MATLAB codes on Intel(R) Core(TM) 17 CPU with 12GB memory and collect the total 
running time of 10^ trials for ID-OMP and 2D-0MP respectively. Because two algorithms output exactly 
the same results, only the speedup of 2D-0MP over ID-OMP with respect to k is reported in Fig. [T] 
From Fig. \T\ we can draw two conclusions: 
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1) As A: increases, the speedup of 2D-0MP over ID-OMP descends gradually. This is because for small 
k, the complexity of 2D-0MP and ID-OMP is dominated by the Project step. As k increases, the 
complexity of other steps will weight heavier. Especially, at the Renew Weights step, the complexity 
of t X t matrix inverse is O(t^), which will ascend quickly as t increases. 

2) As m increases, the speedup of 2D-0MP over ID-OMP becomes more significant. When m = 16, 
the speedup of 2D-0MP over ID-OMP ranges from 10 to 11 times, while when m = 32, the 
speedup of 2D-0MP over ID-OMP ranges from 32 to 35 times. This is because the speedup of 
2D-0MP over ID-OMP comes mainly from the Project step, while at other steps 2D-0MP shows 
Uttle superiority over the ID-OMP. As m increases, the complexity of the Project step weights 
heavier, which explains the aove phenomenon. 

V. Conclusion 

For 2D sparse signal recovery, this paper develops 2D-0MP algorithm. We prove that 2D-OMP is 
equivalent to ID-OMP, but it reduces recovery complexity and memory usage. Hence, 2D-0MP can be 
used as an alternative to ID-OMP in such scenarios as compressive imaging, image compression, etc. 

Following the deduction in this paper, the extension of 2D-0MP to higher dimensional OMP is 
straightforward. For example, by utilizing 3D separable sampling, 3D-0MP can be obtained by defining 
each atom as a 3D matrix. Then at each iteration, the decoder projects 3D sample matrix onto 3D atoms 
to select the best matched atom, and then renews the weights for all the already selected atoms via the 
least squares. 3D-0MP can find its use in hyperspectral image compression. 
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